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A wave equation, that governs finite amplitude acoustic disturbances in a thermoviscous Newto- 
nian fluid, and includes nonlinear terms up to second order, is proposed. In contrast to the model 
known as the Kuznetsov equation, the proposed nonlinear wave equation preserves the Hamilto- 
nian structure of the fundamental fluid dynamical equations in the non-dissipative limit. An exact 
traveling front solution is obtained from a generalized traveling wave assumption. This solution is, 
in an overall sense, equivalent to the Taylor shock solution of the Burgers equation. However, in 
contrast to the Burgers equation, the model equation considered here is capable to describe waves 
propagating in opposite directions. Owing to the Hamiltonian structure of the proposed model 
equation, the front solution is in agreement with the classical Rankine-Hugoniot relations. The ex- 
act front solution propagates at supersonic speed with respect to the fluid ahead of it, and subsonic 
speed with respect to the fluid behind it, similarly to the fluid dynamical shock. Linear stability 
analysis reveals that the front is stable when the acoustic pressure belongs to a critical interval, 
and is otherwise unstable. These results are verified numerically. Studies of head-on colliding fronts 
demonstrate that the front propagation speed changes upon collision. 

PACS numbers: 43.25.Cb, 43.25.Jh, 43.25.Ts 
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I. INTRODUCTION 



The "classical" equation of nonlinear acoustics pQ , the 
so-called Kuznetsov equation [2] , governs finite amplitude 
acoustic disturbances in a Newtonian, homogeneous, vis- 
cous, and heat conducting fluid. The model equation and 
its paraxial approximation, the Khokhlov-Zabolotskaya- 
Kuznetsov (KZK) equation [3J [3J, are occasionally en- 
countered within studies related to nonlinear wave prop- 
agation. See e.g. the recent works by Jordan [3] who 
presented the derivation and analysis of an exact trav- 
eling wave solution to the one-dimensional Kuznetsov 
equation, and by Jing and Cleveland [5] who described a 
three-dimensional numerical code that solves a general- 
ization of the KZK equation, and the references cited in 
the introductory sections of those papers. Other recent 
works based on the Kuznetsov equation include: anal- 
ysis of energy effects accompanying a strong sound dis- 
turbance [B], studies of generation of higher harmonics 
and dissipation based on a 3D finite element formulation 
7 , and studies of nonlinear wave motion in cylindrical 
coordinates |8. . The derivations of the Kuznetsov equa- 
tion U [TU] and related model equations [TTJ [TJJ [T3J 
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are based on the complete system of the equations of 
fluid dynamics. It has been demonstrated that this sys- 
tem of equations is of Hamiltonian structure in the non- 
dissipative limit |14j . However, in the non-dissipative 
limit, the Kuznetsov equation does not retain the Hamil- 
tonian structure. 

In this paper we propose a nonlinear wave equation, 
which, in the non-dissipative limit, preserves the Hamil- 
tonian structure of the fundamental equations. Further- 
more, we present the derivation and analysis of an ex- 
act traveling front solution, which applies equally well to 
the proposed nonlinear wave equation and the Kuznetsov 
equation. The derivation of the exact solution is based 
on a generalized traveling wave assumption, which leads 
to a wider class of exact solutions compared to the one 
reported by Jordan |U [15]. Furthermore, the introduc- 
tion of the generalized assumption is necessary in order to 
interpret the results of numerical simulations of head-on 
colliding fronts presented in this paper. In order to re- 
late our results to the classical literature, we demonstrate 
that the exact front solution retains a number of prop- 
erties of the fluid dynamical shock. The paper is struc- 
tured as follows: The proposed equation and its Hamilto- 
nian structure are discussed in Section|H] Section UHl con- 
tains the derivation of our exact traveling front solution 
and analysis of its stability properties. In Section IV we 



demonstrate that the front is related the classical shock. 
Section |V] presents numerical investigations of the front, 
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while Section IVII contains our conclusions. 



II. NONLINEAR WAVE EQUATIONS 

Equations governing finite amplitude acoustic distur- 
bances in a Newtonian, homogeneous, viscous and heat 
conducting fluid may be derived from the following four 
equations of fluid dynamics: the equation of motion 



P\ ^7 + ( u ' V ) u 



dt 



= -Vp + 7?Au+(^ + C)V(V-u), (1) 



the equation of continuity 

dp 



ot 



V • (pu) - 0, 



the heat transfer equation 



(2) 



TABLE I: Values of Co, B/A, and 6 for three different sub- 
stances. The values for b are rough estimates obtained from 
Eq. Q neglecting the influence of bulk viscosity and thermal 
losses. 



Substance 


c (ms _1 ) 


B/A 


b (m 2 s- L ) 


Water 


1483 (20°C) a 


5.0 (20°C) b 


1.3 x 10~ 6 (20°C) C 


Air 


343 (20°C) a 


0.4 (20°C) b ' d 


21 x 10 -6 (27°C) C 


Soft tissue 


1540 c 


9.6 (37°C) b ' f 


N/A 



a Ref. [19] 
b Ref. OH 

c Values for po and rj are obtained from Ref. 1191 
d Diatomic gas 
c Ref. [201 

f Human breast fat 



where B/A is the fluid nonlinearity parameter [17] and 
Cq = (dp/dp) = is the small-signal sound speed. 
Then, from Eqs. ([1J)-(|6| we obtain the following non- 
linear wave equation 



,^(| + (u-V) S 



duj 



du 3 _ 
dxi 

C(V-u) 2 



2 du k Y 
+ kAT. 



(3) 



and the equation of state 

P = p(p,s)- 



Here x = (x, y, z) are the spatial (Cartesian) coordinates 
and t denotes time, u = (u, q, w) is the fluid particle 
velocity, p is the density of the medium, p, s, and T 
are the thermodynamic variables pressure, entropy and 
temperature, respectively, rj and £ are the coefficients of 
shear and bulk viscosity, and k is the heat conductivity 
coefficient. A is the Laplace operator. 

To obtain a nonlinear wave equation all dependent 
variables except one are eliminated from the system 0- 
Q , resulting in a nonlinear wave equation for that single 
variable M \M HU EE2 US]- The deviations of p, p, 
s, and T from their equilibrium values, po, po, sq, and 
T are assumed to be small, as well as the fluid parti- 
cle velocity, |u|. The heat conductivity coefficient k and 
the viscosities r\ and £ are also treated as small quanti- 
ties. In order to obtain a second order approximation, all 
equations are written retaining terms up to second order 
in the small quantities. It is assumed that the flow is 
rotation free, V x u = 0, thus 



u = -V^, 



(5) 



where tp is the velocity potential. Furthermore, it has 
become customary to use the following approximation 
for the equation of state [16] 



P - Po = Cq (p - Po) + 



cl B/A 



Po 



(p- poY 



dp 
ds 



(s - s ) , (6) 



dip 
~di 




bAiP + {V4>f 




LL 



(4) where b is the diffusivity of sound 



Po 13 \C V 



1 



(7) 



(8) 



and Cy and C p denote the heat capacities at constant 
volume and pressure, respectively. Typical values of the 
physical parameters cq, B/A, and b are given in Table [I] 
In the first order approximation, Eq. Q reduces to 



d 2 ip 



= c 2 Aip. 



(9) 



Introducing Eq. ([9| in the first term on the right hand 
side of Eq. ^ , the Kuznetsov equation [5] 



d 2 %P 



clAiP 



bAiJ; + (Vip) 2 




is obtained. 

In absence of dissipation, i.e. rj = £ = 0, Eqs. and 
([2]) possess Hamiltonian structure [T3] . This property is, 
however, not retained in Eq. ( fTo| ) with 6 = 0, i.e. the non- 
dissipative limit of the Kuznetsov equation is not Hamil- 
tonian. In contrast, Eq. ([7| does retain the Hamiltonian 
structure in the non-dissipative limit. Accordingly, the 
equation may be derived from the Lagrangian density 



C 



B/A -\^f 



64 



(11) 
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using the Euler-Lagrange equation 1 . From the Legendre 
transformation [22] we obtain the corresponding Hamil- 
tonian density as 



H = 4 



2 (W) 2 , (^ B/A-l ^ 



3c 2 , 



2 2 

which may be integrated to yield the total Hamiltonian 

»+oo Z'+OO /*+oo 



H 



— OO J— OO J — OO 



T[ dx dy dz. 



(13) 



and 



P - Po = Po [ipt 



2 



2c, 



\ {A? 



(17) 



respectively. It should be noted that Eqs. (16 1 and (17 1 



are derived from the fundamental equations, thus, the 
expressions are not specific to any of the two model equa- 
tions and ([TO). 



Taking the time derivative of H in Eq. (13 1 with Ti re- 
placed by Eq. (12 1, and using Eq. Q, one can obtain a 
simple expression for dH/dt. Doing this in one spatial 
dimension we obtain after some calculations 



dH 
~dt 



+00 



-b 



(ip xt ) dx. 



(14) 



In Eq. ( 14 1, which is sometimes called the energy balance 



equation, we observe that the first terms on the right 
hand side correspond to energy in- and output at the two 
boundaries, and that the last term accounts for energy 
dissipation inside the system. 

In the remaining portion of this paper we shall limit 
the analysis to one-dimensional plane fields, in which case 
the proposed model equation ^ reduces to 



iptt ~ c tp x 



d_ 
Ft 



btp xx 



B/A-l 



2c 2 



(15) 



where subscripts indicate partial differentiation. 

Finally, for later reference we give the second-order 
expressions for the acoustic density, p — po, and acoustic 
pressure, p — po, in terms of the velocity potential, ip. 
From the equations of motion ([lj and state Q6h , subject 
to the basic assumptions of the derivation of the two 
model equations |7]) and ( |To| ), we obtain 



III. EXACT TRAVELING FRONT SOLUTION 

Recently, a standard traveling wave approach was 
applied to the one-dimensional approximation of the 
Kuznctsov equation ( 10 ) to reveal an exact traveling wave 
solution |4j[15j. In this section we extend the standard 
approach by introducing a generalized traveling wave as- 
sumption and analyze the stability properties of the so- 
lution. 



A. Generalized traveling wave analysis 

We introduce the following generalized traveling wave 
assumption 



ip(x, t) = *(.t - vt) - Xx + at 
= - Xx + at, 



(18) 



where A and a are arbitrary constants, v denotes the 
wave propagation velocity, and £ = x — vt is a wave 
variable. The inclusion of —Xx + at in Eq. (18 1 leads 



to a wider class of exact solutions, compared to the one 
obtained from the assumption ip = ^>(x — vt), which is 
the standard one. Furthermore, the introduction of the 
generalized assumption is necessary in order to interpret 
the results of numerical simulations of head-on colliding 

into 



fronts presented in Section VB 



the nonlinear wave equation ( 15 ) we obtain the ordinary 
differential equation 



Inserting Eq. ( 18 



B/A - 1 a 



tyxxj > (v 2 - cl) 9" = (-vV + a) 9" - v^y>9" 

+ (9'-X) 2 + ^^(-v9' + af }. (!«)) 



1 Letting r) = £ = 0, u = -Vi/ 1 , and p = pT/7 in Eqs. (l]{2) 
respectively, one can derive the Lagrangian density 

7 



CpEE 



and 



corresponding to the potential Euler equation (PEE) given in 
Ref. 1211 Expan ding Cpee to third order and letting 7 = B/A+l 
we obtain Eq. 



where prime denotes ordinary differentiation with respect 
to £. Integrating once and introducing $ = —9', Eq. ( 19 1 
reduces to 



C : rh¥ ['- + B/A ? \ 2 ) r.L-' 
.2 2c 2 



B/ \ 1 v)v 2 -2\v-(%-a}3>, (20) 
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where C is a constant of integration. Requiring that the 
solution satisfy $' — > as £ — > ±00, and either 



where 9 is an arbitrary constant, lead us to C — and 



B ' A - l e^ 



2c 2 



1 > 2 

1 n a V 



+ 2X )v + 4 + a = 0. (22) 



In order to obtain our traveling wave solution, we sep- 
arate the variables in Eq. (20 1 subject to C = 0, then, 



using Eq. (22 1, we find the solution to be the traveling 
front 



46 



2r 2 



(23) 
(24) 



where £0 is an integration constant, \l\ is the front thick- 
ness, and < <& < 9. Finally, using $ = — ty' and in- 



serting Eq. (23 1 into Eq. (18 1 we obtain (apart from an 



arbitrary constant of integration) 



Xx + at, 
(25) 



which is the exact solution for the velocity potential. 

Traveling tanh solutions, such as the front solu- 
tion (23 1, are often called Taylor shocks. The existence 



of an exact solution of this type to the classical Burgers 
equation is a well known result |23j . However, the Burg- 
ers equation is restricted to wave propagation either to 
the left or to the right. The model equation considered 
in this paper does not suffer from this limitation, as shall 
be illustrated in Section PVBI 

Regarding the exact solution derived above, the phys- 
ical properties of the flow associated with the traveling 
front are obtained from the partial derivatives of Eq. (251, 
which are given by 



-ip x = <f> + A and ip t 



(26) 



According to Eq. ^ the fluid particle velocity is ob- 
tained as it = —ipx, and the first order approximation of 
Eq. (17 1 yields the acoustic pressure as p — po « Poipt- 
The boundary conditions of the front are obtained from 



Eqs. (231 and (26 1 as 



9 + x, e- 

A, £- 
v9 + a, £ 



=Foo 
±00 ' 

-> Too 
-> ±00 



(27a) 
(27b) 



where upper (lower) signs apply for I > (I < 0). Hence, 
the four parameters v, 9, A, and a, that was introduced 
in the derivation of the exact solution, determine the four 
boundary conditions of the front. From these boundary 
conditions we find that 9 and v9 correspond to the heights 
of the jump across the front measured in — ip x and ipt, 
respectively, see Fig. [I] At this point it is appropriate to 
emphasize that, in order for the exact solution to exist, 
the four parameters v, 9, A, and a must satisfy the cubic 
equation (22 1. Furthermore, the allowable values of the 



wave propagation velocity correspond to the real roots of 
this equation. A noticeably property of Eq. (22 1, which 



will prove useful later on, is that the equation is invariant 
under the transformation 



X^9 + X, 



(28) 



Also the boundary conditions (27 1 are invariant, since, 



according to Eq. ( 24 ) , the above transformation leads to 
I -> —I. 



1 




v9+a 




FIG. 1: The exact solution of Eq. (15 1 represents a traveling 



front. In order for the solution to exist, the wave propagation 
velocity, v, the front height, 8, and the two constants, A and 
a, must satisfy Eq. (|22[). The plot shows a front with I > 0. 



In order to investigate the relationship between the 
front height, 9, and the front propagation velocity, v, we 
solve Eq. ( 22 1 with respect to 9 to obtain 



1 - B -^ 2 — 1 a ) v 2 - 2Xv -c\ - n- 



B/A- 



(29) 



2r 2 

For B/A < 1 the curve 9(v) has singularities at 

1/2 



v = v* = ± 



1 -B/A 



(30) 



and for B/A > 1 the curve has a maximum 2 at (v, 
(«max,#max), where u max is obtained as 



'"max — CO 



'3 B/A + y/9 (B/A) 2 + 12(B/A - 1) 
2(B/A-1) 



1/2 



(31) 



2 The critical point (v, 9) = ("Umax, #m£ 
[4] as the solution bifurcation point. 



) was identified by Jordan 



5 



when A = and a = 0. These two characteristic proper- 
ties of the curve are illustrated in Fig. [2] 




FIG. 2: The plots show the relationship between the front 
height, 8, and the front propagation velocity, v, given by 
Eq. (29| with A = 0, a = 0, and B/A = {0, 0.8, 1, 1.1, 1.5, 5} 
(see labels on the plots). The dashed lines indicate the singu- 
larity at v = v s , which is defined in Eq. (30 1, and crosses in- 
dicate the maximum (8, v) = (# max , w max ) defined in Eq. (31 1. 



Finally, it should be emphasized that the general- 
ized traveling wave analysis conducted above also applies 
to the one-dimensional approximation of the Kuznetsov 
equation (10). In this case Eq. (22) is replaced by 3 



B/A 



6v A - 1 - 



B/A 



a )v z + (9 + 2A) v + 4 = 0, (32) 



and Eq. (24 1 by 



4b 



B/A 
zl 



(33) 



Apart from these changes, a generalized traveling wave 
analysis of the Kuznetsov equation is basically identical 



to that of Eq. (|15j). 

The Hamiltonian structure, however, is unique to 
the proposed nonlinear wave equation |7| and its one- 
dimensional approximation Eq. (15). In order to estab- 



lish a relationship between the exact solution, derived in 
this section, and the Hamiltonian structure of the gov- 



erning equation, we insert Eq. (25 1 into Eq. (14 1 and the 
one-dimensional approximations of Eqs. (12 1 and (|13l 



Then, after some calculations, we find that Eq. (I14J re- 
duces to the cubic equation (22). Hence, the exact trav- 



eling front solution of the proposed Hamiltonian model 
equation (15) satisfies the energy balance equation (14 1. 



B. Linear stability analysis 

In order to gain insight into the stability properties 
of the traveling front solution we initially consider the 



3 Eliminating A and a from Eq. ( |'S-'| makes the equation equivalent 
to the previously reported result _4. 



constant solution 

-*!>* = K 



and 



(34) 



which satisfies the nonlinear wave equation (151. The two 



constants K and L are arbitrary. In order to investigate 
the linear stability properties of this solution, we add 
small perturbation terms to the constant values as 



-ip x = K - e\x 



and 



ip t = L + ext, 



(35) 



where \ — x{ x >t) an d £« 1. Then, inserting Eqs. (351 
into Eq. ( |15[ ) and keeping terms up to first order in e we 
obtain the following linear perturbation equation 



A B/A-l\ (2 , 

^1 ^ L J Xtt - (c + L) Xxx -- 

Inserting the single Fourier mode 

X (x,t) = De i ( kx -" *>, 



bXxxt - ZKxxt- 
(36) 

(37) 



where D is the amplitude, k is the wave number, and 
u> is the angular frequency, into Eq. ( 36 1 , we obtain the 
following dispersion relation 



(2Kk - ibk 2 ) oj 

+ (eg + L) k 2 = 0. (38) 



The constant solution (34 1 is asymptotically stable only 
if all solutions of Eq. (36 1 approach zero as t — > oo. 



This is the case when the imaginary part of both roots 
in Eq. (38 1, u>% and UJ2, are negative. It can be shown 



that for B/A > 1 the only requirement in order for the 
imaginary part of both roots to be negative is 



-Cn < L < 



B/A- 1 



(39) 



When B/A < 1 the only requirement for both roots to 
have a negative imaginary part is 



— c n < L < oo. 



(40) 



Hence, the stability properties of the constant solution 
(34 1 are determined exclusively by L, i.e. the constant 
value of ipf Recall that the acoustic pressure is pro- 
portional to ift, thus, the level of the acoustic pressure 
determines the stability properties of the solution. 

In order for the front solution to be stable, it is a neces- 
sary condition that both left and right asymptotic values 
of ipt, given by Eq. (27b I, belong to the interval (39 1 



when B/A > 1, and the interval (40) when B/A < 1. 
In Section |"V A| we shall further investigate this stability 
criterion by means of numerical simulations. 
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IV. FRONT-SHOCK RELATIONSHIP 



Furthermore, changes in products of ip x and ipt are 



Within fluid dynamics, a shock denotes a sharp change 
of the physical quantities. A shock propagates at super- 
sonic speed with respect to the fluid ahead of it, while 
it remains subsonic with respect to the fluid behind it. 
The physical quantities of the flow on each side of the 
shock are connected by the Rankine-Hugoniot relations, 
which are conservation equations for mass, momentum 
and energy. In the following we shall demonstrate that 
the front solution of the proposed Hamiltonian model 
equation (15) retains these properties. 



A. The Rankine-Hugoniot relations 

Using square brackets to denote the change in value of 
any quantity across a shock, e.g. 



[p] = P* 



Pb, 



(41) 



where b denotes the value behind the shock and a de- 
notes the value ahead of it, the Rankine-Hugoniot rela- 
tions may be written as [H] 



mass : 
momentum : 

energy : 



[p(u-v)]=0, 
p + p(u — v) 2 = 0, 



h+ (u- vf /2 



0. 



(42) 
(43) 
(44) 



where v is the shock propagation velocity and h is the 
enthalpy. 

We now replace u, p, p, and h with expressions in terms 
of ip x and ipt , and write all equations retaining terms up 
to second order. Upon setting u = —ip x and substituting 
Eqs. JTol) and (fl7| into Eqs. p2l and (E3l we thus obtain 



B/A-l 2 \ 



2cl 



■ Ipt {i>x +V) - cllpx 



0, (45) 



and 



B/A-l 



(iptf v 2 - ip t v z - 2tp t tpxV 



2cl 



(V* 



2clip x v - clipt 



0, (46) 



respectively. The dissipative terms involving n, £, and 77 
do not appear in Eqs. (45 I and (46 1, since ip xx — * ahead 



of and behind the front. The changes in ip x and ipt across 
the front are obtained from the boundary conditions ( 27 1 . 



Assuming that I > and v > 0, and using the notation 
introduced in Eq. (41 1 we may write 



(47a) 



Ipxlpt 



= -9 2 - 26X, 
= -v 2 9 2 - 2v6a, 
= v8 2 + v8\ + da. 



(47b) 
(47c) 
(47d) 



Inserting Eqs. (47 1 into Eqs. (45 1 and (46 1 both conser- 
vation equations reduce to the cubic equation ( 22 1 . This 



striking result leads to the conclusion, that Eq. (22 1 im- 
plies conservation of mass and momentum. At this point 
it should be noted that the generalized traveling wave 



analysis of the Kuznetsov equation ( 10 1 leads to the cu- 



bic equation (32 1, which is not in agreement with the 



conservation equations for mass and momentum. 

In order to handle the enthalpy in the condition for en- 
ergy conservation ( |44| we shall make use of the following 
fundamental thermodynamic relationship |24j 



Vh 



Vp 



(48) 



Using the equation of motion Q, subject to the basic 
assumptions of the derivation of the model equations in 
Section [Tl] we obtain from Eq. (j44j) 

[ipt + vip x ] = 0. (49) 

Alternatively, Eq. p9"] ) follows directly from the general- 
ized traveling wave assumption ( 18 1. Hence, the traveling 



wave assumption implies energy conservation in the flow. 



B. Sub-/supersonic speeds of propagation 

In order to determine whether the traveling front so- 
lution, derived in Section [ill A[ propagates at sub- or su- 
personic speed with respect to the fluid ahead of it and 
the fluid behind it, we need to introduce the speed of 
sound in these regions of the fluid. Without loss of gen- 
erality, we may consider only fronts propagating in the 
positive direction, v > 0, since Eq. ( 15 ) is invariant under 
the transformation x —x. Furthermore, we shall limit 
the analysis to stable fronts, i.e. ip t must belong to the 



interval (39) when B/A > 1, and the interval (40) when 
B/A < 1. Then, letting 9 -> in Eq. p2j and solving 



for v yields the small signal propagation velocity, which 
is equivalent to the speed of sound, c. Introducing A = K 
and a = L we obtain 



V = c{K, L) 




K 2 + (L + c§) I 1 



B/A-l 



L 



(50) 



where K and L denote the constant levels of — ip x and tp t 



respectively, at which the speed of sound ( 50 ) is evalu 



ated. Inserting the boundary conditions of the front into 



7 



Eq. (50 1, i.e. substituting Eq. (27a I for K and Eq. (27b) Eq. (57). Non-dimensional versions of the parameters, v, 



for L, we obtain the speed of sound ahead of, c a , and 
behind, Cb, the front 



c a = c(A,ct), 

b 

c b = c(9 + X,v 8 + a), 



(51) 
(52) 



where upper (lower) subscripts apply for I > (I < 0). 
Note that, under the transformation (28), Eq. (52 1 trans- 



forms into Eq. (51 1. Hence, without loss of generality we 



shall consider only Eq. ( |51[ ) in the following. 

In order to compare the front propagation velocity, v, 
to c a and Cb we make the following observations. Insert- 
ing Eq. d29} into Eq. Q yields 



I 



Abv 



B/A-l 



(53) 



The denominator in Eq. ( 53 ) becomes zero when v = 



c(X,a), where c(X,a) is given by Eq. (50 1. Then, given 
that v > 0, we obtain from Eq. (53 1 that 



v > c(A, a) <^ I > and v < c(A, a) ^ I < 0. (54) 



Finally, from Eqs. (51 1 and (54 1 it follows that 
v > c a and v < c^. 



(55) 



Hence, in all cases, the propagation velocity of the exact 
traveling front solution is supersonic with respect to the 
fluid ahead of the front, and subsonic with respect to the 
fluid behind it. 



V. NUMERICAL RESULTS 

All numerical calculations rely on a commercially avail- 
able software package 4 , which is based on the finite ele- 
ment method. For convenience we introduce the follow- 
ing non-dimensional variables, denoted by tilde 



i==±t. (56) 



Under this transformation we may write Eq. (15) as 



ipu - i^xx = iptip. 
d 



dt 



4>xx + (ipxY 



(57) 



where the tildes have been omitted. From a comparison 
of Eqs. (151 and (57 1, we find that the results of the 
previous sections subject to b = 1 and Cq = 1 apply to 



4 COMSOL version 3.2a, http://www.comsol.com (2005) 



A, and a, also indicated by tilde, become 
A a ~ 9 



A 



Co 



Co 



CO 



(58) 



In the following analysis we consider only the non- 
dimensional formulation of the problem. For notational 
simplicity we shall omit the tildes. 



A. Investigation of the front stability criterion 

In order to investigate, numerically, the stability prop- 
erties of the front, we chose as initial condition for the 



numerical solution, the exact solution given by Eqs. (25 1 
and (26 1, and choose v, 9, A, a, and B/A such that 
Eq. (22) is satisfied. For the sake of clarity, we shall 



limit the numerical investigations to the specific case of 
A = 0, cr = 0, v > 0, and I > 0, which, according to 



Eq. (27 1, corresponds to fronts propagating to the right 



into an unperturbed fluid. 

A first numerical simulation is presented in Fig. [3] Ev- 
idently, the numerical algorithm successfully integrates 
the initial condition forward in time. This finding indi- 
cates that, for the specific choice of parameters, v = 1.4, 
9 = 0.127, and B/A = 5, the front exists and is stable. A 
second numerical simulation is presented in Fig. [4] This 
initial condition is given a larger velocity, v = 1.7, and 
a larger height, 9 = 0.153, while B/A = 5 remains un- 
changed compared to the first example. The parameters 
are chosen such that Eq. (22 1 remains satisfied. Appar- 



ently, the numerical algorithm fails when integrating the 
solution forward in time, which indicates that the front 
is unstable for the specific choice of parameters. Given 
that a = and I > 0, the left and right asymptotes of 
the front are given by ip t — v9 and ipt — 0, respectively, 
according to Eq. (27b). Clearly, the right value belongs 
to the interval (39 1, thus, it does not causes instability 
of the front. However, if the left value, v8, lies outside 



the interval (39), it causes instability of the front. Insert- 



ing the value of B/A from the two examples above into 
Eq. (39), we find that in the first and second example, 
v9 lies inside and outside the interval ( |39| , respectively. 
Hence, the behavior observed in Figs . [3] and [4] agre es with 
the stability criterion introduced in Section |III B| 

A large number of numerical simulations have been 
performed in order to systematically investigate the sta- 
bility properties of the front. Within each simulation 
the parameters in the initial condition are, again, chosen 
such that Eq. (22 1 is satisfied. The result of this inves- 
tigation is presented in Fig. [5] Still, a = and I > 0, 
such that the left asymptotic value of the front is given 
by ipt — vO. For B/A > 1 the stability threshold curve in 
the (B/A, ii)-plane is obtained when v9 equals the upper 
bound of the interval (|39|. Using Eq. (p9| we obtain 



1 



B/A-l 



y/(B/A-l)(2B/A+l) 
^ V= BjA—l ' (59) 



8 



For B/A < 1, the stability threshold curve is given by 
Eq. (30), since v6 lies within (outside) the interval (40) 
when v < v s (v > t> s ), according to Eq. d29b. The two 
stability threshold curves are included in Fig. 151 The fine 
agreement between the numerical results and the stabil- 
ity threshold curves indicates that the stability criterion, 
introduced in Section |III B| is both necessary and suffi- 
cient in order for the front solution to be stable. 




FIG. 3: The initial condition at t = (bold line) is obtained 
from Eqs. Q and (26]} subject to v = 1.4, = 0.127, A = 0, 
a = 0, and B/A = 5, which satisfy Eq. (22 I. The numerical 



solutions are shown over the time interval < t < 25. 
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FIG. 4: The numerical algorithm fails to integrate this solu- 
tion forward in time (insert shows magnification). See caption 
of Fig. [3] v = 1.7, 6 = 0.153, A = 0, a = 0, and B/A = 5. 
The numerical solution is shown at time t = 3.8 x 10 -3 . 
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FIG. 5: Each point in the (B/A, v) -plane represents a nu- 
merical simulation, with the initial condition obtained from 
Eqs. (25 I and (261 subject to A = 0, a = 0, and 9 given by 
Eq. (29 1. Crosses and circles indicate stable and unstable so- 
lutions, respectively (compare with Figs. [3] and [4|. Solid lines 
represent the stability threshold curves given by Eqs. ( |30[ ) and 
(591. 



B. Head-on colliding fronts 

The numerical simulation presented in Fig. [6]shows the 
result of a head-on collision between two fronts. From the 
simulation we observe that two new fronts emerge upon 
the collision. The contour plot reveals that these fronts 
travel at a higher speed, compared to the speed of the 
fronts before the collision. For other choices of initial 
condition, we found the outcome of the head-on collision 
to be fronts traveling at lower speed, compared to that 
of the fronts before the collision. 

In order to analyze solutions of Eq. ( 15 1 that com- 



prise two fronts, we assume that these fronts belong to 
the class of exact front solutions derived in Section IIII Al 
above. Investigations of the fronts that emerge upon a 
head-on collision have made it clear that this assumption 
is true, only when the generalized traveling wave assump- 
tion is considered, in contrast to the standard traveling 
wave assumption. Then, for each of the two fronts in the 
solution we introduce four new parameters, v, 0, A, and 



a, which must satisfy Eq. (22 1 as 



B/A- I 



-(l-(B/A-l)ai)vf 
9i + 2\i J «i + <Ti + 1 = 0, i = l,2, (60) 



where subscript 1 and 2 denote parameters associated 
with waves positioned to the left and right, respectively. 
Furthermore, we require that solutions comprising two 
fronts are continuous and satisfy the following set of ar- 
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bitrary boundary conditions 

/ P, x — > -00 , 



i?, x — * — oo 
S, a; — > +00 



(61) 



Assuming that l\ > and Z2 < 0, we find, using Eq. (27 1, 



that the these requirements lead to the following condi- 
tions 

\i = \2 = P-0 1 = Q- 9 2 , (62a) 
a 1 = <T 2 = R-v 1 6 1 = S-v 2 6 2 , (62b) 
61 + 6 2 = P - Q, v x 6i + v 2 6 2 = R~S. (62c) 

Then, we substitute the boundary values found in Fig. [6] 
for P, Q, R, and S in Eqs. (62 1, and substitute the 



value of B/A into Eq. (60 1. Finally, solving the system 
of equations p0| and (62 1, we obtain the results listed 
in Table [TTJ The solution in the first row of the table 
corresponds to the initial fronts found in Fig. [6] The 
solution in the second row corresponds to two unstable 
fronts, according the stability criterion discussed above. 
The two fronts that emerge upon the head-on collision 
are defined by the values found in the third row of the 
table. Hence, the fronts after the collision travel at the 
velocities —v\ = v 2 = 1.76, which is in agreement with 
the velocities determined from the slope of the contour 
lines in Fig. [6] 



VI. CONCLUSIONS 

A nonlinear wave equation that governs finite ampli- 
tude acoustic disturbances in a thermoviscous Newtonian 
fluid, and includes nonlinear terms up to second order, 
has been presented. The single dependent variable is 
the velocity potential. It has been demonstrated that, 
in the non-dissipative limit, the equation preserves the 
Hamiltonian structure of the fundamental fluid dynam- 
ical equations, hence, the model equation is associated 
with corresponding Lagrangian and Hamiltonian densi- 
ties. Furthermore, we found that the Kuznetsov equation 
is an approximation of the proposed nonlinear wave equa- 
tion. However, in the non-dissipative limit the Kuznetsov 
equation is not Hamiltonian. Exact traveling front solu- 
tions, for the partial derivatives with respect to space and 
time of the dependent variable, has been obtained using 
a generalized traveling wave assumption. This general- 
ized assumption leads to a wider class of exact solutions 
compared the one obtained from a standard traveling 
wave assumption, since the generalized assumption in- 
cludes two arbitrary constants, which are added to the 
partial derivatives. As a result of the generalized trav- 
eling wave analysis we found that, in order for the front 
to exist, its boundary values, its propagation velocity, 
and the physical parameters of the problem must satisfy 
a given cubic equation in the front propagation veloc- 
ity. The derivation of the exact solution applies equally 




100 




FIG. 6: The initial condition (bold lines in the two topmost 
plots), corresponds to two fronts that make a head-on collision 
at t — 42. The initial fronts are defined by v\ — —V2 = 1.19, 
0i = -02 = 8.07 x 10~ 2 , Ai = A 2 = 0, 0-1 = cr 2 = 0, and 
B/A — 5, where subscript 1 and 2 relate to fronts positioned 
the left and right, respectively. For each of the two fronts the 
parameters satisfy Eq. ( |22[ |. Lowermost: contour lines given 
by —ipx = Z, where Z takes 4 equidistantly spaced values 
across each front. 



well to the proposed Hamiltonian model equation and 
the Kuznetsov equation. Results for both equations have 
been given. 

It has been demonstrated that the overall stability 
properties of the front are determined by the stability of 
the two asymptotic tails of the front. A linear stability 
analysis of these steady parts of the solution revealed that 
the front is stable when the partial derivative with respect 
to time, which is proportional to the acoustic pressure, 
belongs to a critical interval, and is otherwise unstable. 
This stability criterion has been verified numerically, by 
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TABLE II: 


Solution of Eqs. ( 60 1 and ( 62 1 subject to P 


= -Q = 


8.07 x 10~ 2 , R = 


S = 


9.60 x 10 


" 2 , and B/A = 


5 (compare 


with Fie- Ifi 
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Solution 


V! 6>i 


V2 


e 2 




Ai 


= A 2 


01 = 02 


1 


1.19 8.07xl0 -2 


-1.19 


-8.07x10" 


■i 










2 


-3.25 8.07xlCT 2 


3.25 


-8.07x10" 


2 







35.8xl0~ 2 


3 


-1.76 8.07xl0~ 2 


1.76 


-8.07x10" 


2 







23.8xl0~ 2 



using the exact front solution as initial condition in a 
number of numerical simulations. 

It has been demonstrated that, in all cases, the front 
propagates at supersonic speed with respect to the fluid 
ahead of it, while it remains subsonic with respect to 
the fluid behind it. The same properties have been re- 
ported for the classical fluid dynamical shock. Further- 
more, it has been demonstrated that the cubic equation, 
mentioned above, is equivalent to the well established 
Rankine-Hugoniot relations, which connect the physical 
quantities on each side of a shock. However, this result 
was accomplished only when considering the cubic equa- 
tion obtained from the analysis of the proposed Hamil- 
tonian wave equation. The generalized traveling wave 
analysis based on the Kuznetsov equation is not in agree- 
ment with the Rankine-Hugoniot relations. Estimates of 
the front thickness may be obtained using the values for 
the diffusivity of sound listed in Table. [I] In water and 
air front thicknesses are found to be of the order 10~ 9 
and 10~ 7 meters, respectively However, caution should 
be taken with these estimates, as the small length scales 
violates the continuum assumption of the governing equa- 
tions. 

Numerical simulations of two head-on colliding fronts 



have demonstrated that two new fronts emerge upon the 
collision, and that these fronts, in the general case, travel 
at speeds, which are different from the speeds of the 
fronts before the collision. It has been demonstrated that 
the velocities of the fronts after the collision may be cal- 
culated, based on information about the fronts before the 
collision. However, in order to accomplish this calcula- 
tion, it has proven necessary to introduce the generalized 
traveling wave assumption in the derivation of the front 
solution. 

In future studies, it would be rewarding to further in- 
vestigate a variety of interacting fronts, other than the 
head-on collision reported in this paper. Also a search 
for other types of wave solutions, might learn us more 
about the properties of the proposed Hamiltonian model 
equation and the Kuznetsov equation. 
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